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Abstract 

The periodic, inverse scattering transform (PIST) is a powerful analytical tool in the theory of integrable, nonlinear evolution 
equations. Osborne pioneered the use of the PIST in the analysis of data form inherently nonlinear physical processes. In particular, 
Osborne's so-called nonlinear Fourier analysis has been successfully used in the study of waves whose dynamics are (to a good 
approximation) governed by the Korteweg-de Vries equation. In this paper, the mathematical details and a new application of the 
PIST are discussed. The numerical aspects of and difficulties in obtaining the nonlinear Fourier (i.e., PIST) spectrum of a physical 
data set are also addressed. In particular, an improved bracketing of the "spectral eigenvalues" (i.e., the ±1 crossings of the Floquet 
discriminant) and a new root-finding algorithm for computing the latter are proposed. Finally, it is shown how the PIST can be used 
to gain insightful information about the phenomenon of soliton-induced acoustic resonances, by computing the nonlinear Fourier 
spectrum of a data set from a simulation of internal solitary wave generation and propagation in the Yellow Sea. 
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1. Introduction 

The Korteweg-de Vries (KdV) equation is arguably the most famous "soliton-bearing" equation, governing phe- 
nomena as seemingly disparate as the motion of lattices, the collective behavior of plasmas and the evolution of 
hydrodynamic waves (see, e.g., [1 2 18 20 23 26 28 1 and the references therein). From a mathematical point of view, 
the most interesting property of the KdV equation, discovered by the celebrated numerical experiment of Zabusky 
and Kruskal (ZK) l26l , is the ability of its (localized) traveling-wave solutions to retain their "identity" (shape, speed) 
after colliding with each other. This has been a topic subject to much research over the last few decades, giving birth 
to the idea of integrable equations and leading to the discovery of the inverse scattering transform. In particular, it has 
been shown that the KdV equation is "fully-integrable" on both the real line and periodic intervals; therefore, it can 
be solved exactly by the (periodic) inverse scattering transform [ 1 1 . What is more, Salupere et al. [ 23 ] have shown that 
the ZK experiment features more than just the particle -like interactions of the localized solutions of the KdV equation: 
there exist hidden solitons, and soliton ensembles emerge on long time scales in the "wave soup" created by the ZK 
sine wave initial condition. 
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It is this richness of the solution of the periodic KdV equation and its integrability by the scattering transform that 
form the basis of Osborne's nonlinear Fourier analysis [12|. The latter has been successfully employed in the Fourier- 
like decomposition of data from inherently nonlinear physical phenomena such as shallow-water ocean surface waves 
E2l . laboratory -generated surface waves [20 1 and internal gravity waves in a stratified fluid [28]. In the present work, 
we focus on the implementation and application of Osborne's nonlinear Fourier analysis to yet another nonlinear-wave 
phenomenon governed by the KdV equation: internal solitary waves in the ocean. 

In their 1980 article, Osborne and Burch [18| put forth a model of internal solitary waves in the ocean based on 
the KdV equation. Since then, their model has been applied to a variety of internal solitary wave phenomena (see, 
e.g., 12) and the references therein). In particular, the latter model was employed by Zhou et al. [27] to explain the 
anomalous acoustic signal loss in the presence of internal solitary waves in the ocean. More recently, Chin-Bing et 
al. (6 1 and Warn-Varnas et al. M24I25I have provided evidence of such anomalous signal losses via coupled ocean- 
acoustic simulations. However, though there is ample evidence of the effects of internal solitary waves on acoustic 
signals, there remain unanswered (theoretical) questions regarding, e.g., energy transfer (mode conversions) in the 
acoustic field (see, e.g., [24] and the references therein). Thus, it is a goal of the present work to lay the foundation for 
the use of the nonlinear Fourier analysis of internal solitary wave trains in the ocean in obtaining more accurate and 
"natural" wavenumber (or mode) ranges than the ordinary Fourier transform. These, in turn, can be used to prove or 
disprove conjectures regarding energy transfer in the acoustic field. 

This paper is organized as follows. In Section [2] the KdV-based mathematical model of internal solitary waves in 
the ocean is described. In Section [3] the scattering transform for the KdV equation is presented and its interpretation 
as a nonlinear generalization of the Fourier transform is discussed. In Section [4] the numerical implementation of the 
PIST is addressed, some of the shortcomings of previous algorithms are pointed out, and improvements are suggested. 
In Section|5] the PIST is applied to the analysis of simulated internal solitary wave data and the results are presented. 
Finally, in Section|6| we discuss the future of the present technique. 

2. A mathematical model of internal solitary waves in the ocean 

In this section, we summarize Osborne and Burch's [18| model of internal solitary waves (or solitons, as the case 
might be) in the ocean. It is based upon the assumption of a two-layer stratification in the ocean, where the two 
layers are separated either by the thermocline or by the pycnocline, i.e., a region of rapid change in the temperature 
or the density of the stratification, respectively. Then, it is supposed that solitary waves traveling on the interface 
between the two layers are governed by the ubiquitous KdV equation, whose solitary wave solutions are solitons (so, 
the terms 'solitary wave' and 'soliton' become interchangeable henceforth). The (displaced) interface is a surface of 
constant temperature or constant (potential) density, and it is thus referred to either as an isotherm or as an isopycnal, 
depending on the type of stratification. This physical set up is illustrated in Fig.[T](see also Fig. 3 in |[T8l ). Moreover, 
for the remainder of this paper, we restrict ourselves to the case of a desnity-stratified ocean, so we use the terms 
'pycnocline' and 'isopycnal' instead of 'thermocline' and 'isotherm.' 

Under the assumption that the internal solitary waves are "long" and the top layer of the stratification is "shallow" 
(this will be made more precise below), the governing (KdV) equation of the isopycnal displacement, which we denote 
by rj(x,t), is 

rit + c ri., + ariri x + l3ri xxx = 0, (x,t) e [0,L] x (0,°°), (1) 

where L(> 0) is the length of the domain, and the subscripts denote partial differentiation with respect to an inde- 
pendent variable. In addition, co(> 0), a(< 0) and j3(> 0) are (constant) physical parameters (see, e.g., [2| for their 
interpretation). They can be written in terms of the variables in Fig.[T]as follows ifTSI : 

/ f pi-pi\ ( hih 2 \ -3c (h 2 -hi\ 1 

where h\ < I12 and p\ < p2- Furthermore, we are interested in the periodic Cauchy problem (i.e., the initial-value 
problem subject to periodic boundary conditions). Thus, we have that 

(ri( Xl 0) =tj (x), xe[0,L], 
\r]{x+L,t) = ri{x,t), (*,r)e[0,L]x[0,oo). 
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Fig. 1. Schematic of the structure of internal solitary waves in the ocean. The notation is as follows: pi and pi denote the fluid densities in the 
top and bottom layers, respectively; h\ and h% denote the distances from the unperturbed isopycnal to the free surface and to the ocean bottom, 
respectively; c is a characteristic propagation speed of the solitary-wave packets in the direction indicated by the arrow. 

In practice, however, data is rarely periodic, so considering quasiperiodic (or, more precisely, "almost periodic") 
solutions, i.e., those for which V5 > 3£ = £(5) such that \T](x + £,t) -T](x,t)\ < 8 for all (x,t) G [0,L] x [0,°°), 
might be more appropriate. Though the details of the PIST carry over to this case, difficulties emerge in the physical 
interpretation of the results lfl4ll . Therefore, we restrict ourselves to periodic data sets and solutions. 

Here, we note that more sophisticated continuous stratification models (still in conjunction with the KdV equa- 
tion) have been proposed (see, e.g., (2 | and the references therein) and successfully used in practical modeling [24- 1 . 
However, for the purposes of this work, the well-tested model of Osborne and Burch is sufficient. 

Finally, it is important to outline the limits of applicability of the KdV equation to stratified fluids. The assumption 
of "long, shallow-water" waves we made above can be made more precise by requiring that (see B2I18I1 ) the largest 
wave amplitude is much smaller than the top layer's depth, i.e., [max,-., r)(x,t)]/hi <C 1, and that the characteristic 
width of the waves is much greater than the top layer's depth, i.e., h\/W <C 1, where W can be taken to be, e.g., the 
largest half-width of the waves. In practice, as illustrated by the results below, the former condition is far more difficult 
to satisfy than the latter. 

3. Interpretation of the PIST as a "nonlinear Fourier transform" 

Next, we turn to the relationship between the periodic, inverse scattering transform (PIST) and the (ordinary) Fourier 
transform, and the interpretation of the former as a nonlinear generalization of the latter. To this end, first we note that 
the solution strategy by the scattering transform can be split into two distinct steps: the direct problem and the inverse 
problem. The former consists of solving the eigenvalue problem 

d 2 

Ji"W = E\if, 3>?:=-j^-kr\ {x), (4) 

where X := a/(6/3) is a nonlinearity-to-dispersion ratio, and E G M is a spectral eigenvalue such that \[E g C is a 
wavenumber. Quite serendipitously, Eq. Q has been studied extensively in the literature: in quantum mechanics, 
is the celebrated (time-independent) Schrodinger operator, and, in the theory of ODEs, it is known as Hill's operator (a 
special case of the well-known Sturm-Liouville operator) [5 |. For periodic "potentials," i.e., when T70 (x) = r)o(x + L) 
Vx G [0, L] as we have assumed, it is well-known that the spectrum of the operator is divided into two distinct sets 
depending upon the boundary conditions imposed on the eigenfunctions \j/ 1151111 . Thus, it is common to classify the 
spectral eigenvalues (often referred to as the "scattering data") as belonging to either the main spectrum, which we 
write as the set {£ , j} 2 j=f l , or the auxiliary spectrum, which we write as the set {/i'?}^ =1 , where N is the number of 
degrees of freedom (i.e., oscillations modes). 
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On the other hand, the inverse problem consists of constructing the nonlinear Fourier series from the spectrum 
U {l^j} using Abelian hyperelliptic functions [12 1 1 1411 or the Riemann ©-function 141161 . In former case, which is 
the so-called jj. -representation of the PIST, the exact solution of Eq. ([TJ, subject to the initial and boundary conditions 
given in Eq. ([3]), takes the form 

1 ( N 2N+1 ") 

It is important to note that all nonlinear waves and their nonlinear interactions are accounted for in this linear 
superposition. Unfortunately, the computation of the nonlinear oscillation modes (i.e., the hyperelliptic functions 
{jj,j(x,t)}^ =l ) is highly nontrivial; however, numerical approaches have been developed ll2D and successfully used 
in practice [22 20 1. In addition, we note that the auxiliary spectrum, often referred to as the hyperelliptic function 
"phases," is such that ^ = £1/(0,0) 1151201 . 

Several special cases of Eq. <|5j offer insight into why the latter is analogous to the ordinary Fourier series. In 
the small-amplitude limit, i.e., when max 1)( <C 1, we have fj,,(x,t) ~ cos(kjX — C0jt + <j>j), where kj is the 

wavenumber, (Oj is the frequency and 0; the phase of the mode. Therefore, if we suppose that all the oscillation modes 
fall in the small-amplitude limit, then Eq. Q reduces to the ordinary Fourier series! This relationship is more than 
just an analogy, a rigorous derivation of the (ordinary) Fourier transform from the scattering transform, in the small- 
amplitude limit, is given in ff3l . Next, if there are no interactions, e.g., the spectrum consists of a single wave (i.e., 
N = 1), we have ji\{x,t) = c\^{k\x — a>\t + (j>\ \m\), which is a Jacobian elliptic function with modulus m\. In fact, it 
is the well-known cnoidal wave solution of the (periodic) KdV equation JTJ . 

For the hyperelliptic representation of the nonlinear Fourier series, given by Eq. d5J, the wavenumbers are com- 
mensurable with those of the ordinary Fourier series, i.e., kj = 2nj/L{\ < j < N) [ 14 1. However, this is not the only 
way to classify the nonlinear oscillations. One can use the modulus mj, termed the "soliton index," of each of the 
hyperelliptic functions, which can be computed from the discrete spectrum as 

"' - ^+1-^ (6) 



<%2j+\ — <%2j- 

Then, each nonlinear oscillation mode can be placed into one of two distinct categories based on its soliton index: 

(i) mj > 0.99 => solitons, in particular, cn 2 (x\m = 1) = sech 2 (x); 

(ii) mj < 1.0=^ radiation, in particular, cn 2 (x|m = 0) = cos 2 (x). 

For mj not in either distinguished limit above, the qualitative structure of the nonlinear oscillations is not immediately 
obvious. Boyd [3 | argues that for moderate moduli the polycnoidal wave solutions of the KdV equation are actually 
well-approximated by the first few terms in their Fourier series, showing they are, in fact, not much different from the 
linear waves of the mj < 1.0 limit. 

Furthermore, it can be shown [22 20 1 that the amplitudes of the hyperelliptic functions are given by 



(£ re f — $ij ) , for solitons ; 



{$lj+i — <*2j)> otherwise (radiation) 



where E m f = &ij*+i is the soliton reference level with j* being the largest j for which m ; > 0.99. Then, clearly, the 
number of solitons in the spectrum is A^ so i = j*. 



4. Numerical implementation of the PIST 

In the present work, we have implemented a modified version of Osborne's automatic algorithm [ 15 1 for computing 
the spectral eigenvalues, which we describe in this section. To this end, we begin the solution of the periodic-potential 
eigenvalue problem, given by Eq. by introducing an "appropriate" basis of eigenfunctions and appealing to Flo- 
quet's theorem [5 1 to write the eigenfunctions one period ahead as a linear combination of the current ones: 

$(x+L,E) = a(x+L,E)<f>(x,E) Vxe[0,L], (8) 
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where a(x + L,E) is the monodromy matrix, and 4> is the fundamental matrix of linearly-idependent eigenfunctions. 
That is to say, a "advances" the eigenfunctions one period to the right, and 3* = I . for some eigenfunction 0, 

\o- 0;) 

where a star superscript denotes complex conjugation. Then, the spectral eigenvalues can be related to the monodromy 
matrix (see, e.g., II14I15I and the references therein) as follows: 

Main Spectrum: Sj such that \ tra(x, Sj) = ±1 , (9) 
Auxiliary Spectrum: ju° such that 0C2i(x,Hj) = 0. (10) 

Now, this may seem like a circular definition because we need to know the eigenfunctions (and eigenvalues) to find 
the monodromy matrix. However, it is possible to find a numerically-accessible quantity that we can use to determine 
a without explicitly finding the eigenfunctions. Then, we can compute the spectrum {Sj} U {[ij} of our data set, i.e., 
of the potential — Xt]q(x). 

To this end, we rewrite the Schrodinger eigenvalue problem, given by Eq. Q, as a first-order system: 

d x V(x,E)=B(x,E) x V(x,E), B(x,E)=\ I, (11) 




dx 

where \P := (y/, \j/ x ) T and q(x,E) := Xrio(x) + E. Integrating Eq. ( fTTj ) with respect to x (while keeping E fixed and 
treating it as a parameter), we have 



¥(x,E) = ex Py B(x,E)dxj, (12) 

where xq S [0,L] is an arbitrary base point lfT5l ; here, we set xq = without loss of generality. 

Now, T7o(jc) is usually a discrete data set, i.e., its values are only known for x =X{'.= iAx, where i is an integer such 
that < i < M— 1, M is the number of points in the data set and Ax:= L/M. Therefore, for such a piecewise-constant 
potential [ 13 1, we immediately have [from Eq. ( fT2) )] that 

y(x i+u E) = e A -* B ^- E ^(x i ,E) Vi. (13) 

Iterating the latter relation gives ^(xq+L^E) — M (xq + L,E) 1 ¥(xq,E), where 



M(x + L,E):= Yl e^ 3 ^-^ (14) 

i=M- 1 

is the so-called scattering matrix. Also, the matrix exponential above can easily be computed: 



e A r B(*„£) = 



cos^AxV 'q{xi,E)\ sinfA^V 'q{xi,E)^ j \/q(xt,l 
^-y / q(x i ,E)sm^Axy / q(x i ,E)j cos^AxV^i7E)) 



Vi. (15) 



Furthermore, notice that y / q(xi,E) is either real or purely imaginary, meaning that the above matrix is always real. 
This observation allows for the implementation of the entire algorithm using real arithmetic, speeding it up by a factor 
of four ifTBI . Finally, it can be shown that jtra = jtriVf and Cfei = M21 1151201 . Thus, the scattering matrix M 
is a numerically-accessible analogue of the monodromy matrix a that does not require a priori knowledge of the 
eigenfunctions of J4f. So, we can use M to compute the spectral eigenvalues. 

Unfortunately, even with this knowledge, computing the spectrum {<£}}U {/^}, i-e-, finding the ±1 crossings of 
the Floquet discriminant A(E) := i tr M(xq +L,E) and the zero crossings of M2\{E), is no easy task. Osborne ifTSl 
suggests using an "accounting function" such as 

Af-2 

SUE)] = I I |sgng -sgn[| (16) 

1=0 

where sgn(j) = %/\x\ f° r a ll rea l 1^0 and it vanishes otherwise, to count the zeros of some function ^(E)[— 
M 12(E), A(E), etc] at a certain energy level E. Unfortunately, the accounting function S[ ■ ] suffers numerical instability 
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for the spectral scales that result from the internal solitons model (notice the spike in 5[A(£)] in Fig. [5] for E £ 
[—0.00045,-0,0004]). The latter instability makes the use of such accounting functions undesirable when A 1. 
For example, in the study of internal solitary waves in far deeper water columns than considered here, the nonlinearity- 
to-dispersion ratio is |A| ~ 10~ 8 iflOl . and 5[A(£)] is dominated by numerical noise. 

Therefore, in the present work, we take a more pragmatic approach to computing the spectral eigenvalues. Rather 
than using the automatic algorithm, which (in theory) assures that all eigenvalues are found, we simply define a fine 
enough grid a priori, e.g., {E m \ n + kAE} k ^ Q with AE := (E max — E m i n ) /Ne, where Ne is the number of intervals and 
Emax, £min are appropriately chosen upper and lower bounds on the spectral eigenvalues (see below), so that all zero 
crossings of some function | (E)[— Mn(E), A(E), etc] are resolved. Then, by sampling %(E) on the latter grid and 
checking for sign changes, we obtain a bracketing of its roots. 

Now, the automatic algorithm given in ifTBI requires using S[M\ i(E)] to bracket the zeros of M\ \ (E), which along 
with the zeros of M2i(E) give a bracketing of the ±1 crossings of A(E). However, since the number of zeros of M21 (E) 
is equal to the number of zeros of M\ \ (E) , which is, in turn, equal to the number of oscillations modes in the nonlinear 
Fourier series [ 15], we only compute the zeros of M21 (E), thereby reducing the computational cost of the algorithm. 
Also, at this point, we recall the following important facts stemming from the structure of the equations [ 1 ]: 

(i) M21 (E) = exactly once for E £ [a,b], where a and b are such that A(a) = A(b) = +1 or A(a) = A(b) = — 1, 

(ii) A(E) = exactly once for E £ [a,b], where a and b are such that A(a) = +1 and A(b) = — 1 or A(a) = — 1 and 
A(6) = +l. 

This means that once we have computed all the zeros of M2\(E) and A(E) in [£mm,£max], and ordered them in 
ascending order, we have a complete bracketing of the ±1 crossings of the Floquet discriminant. 

Given the highly oscillatory nature of the A(E) and M2\(E) (see Fig. [5]), keeping the roots bracketed during the 
computations is critical to the success of the algorithm. Originally, Newton's method was used in [20 1, but, due to the 
latter's tendency to "jump out" of the initial bracket, the bisection method was employed in |Q3). The latter guarantees 
that the new approximation to the root remains in the initial bracket, but converges only linearly, unlike Newton's 
method. To get the "best of both worlds," in the present work, we employ the (modified or Illinois-type) regula falsi 
method [7 1. The latter guarantees that successive approximations of the root remain in the initial bracket and converges 
superlinearly. Now, let us summarize our implementation of the regula falsi method for the reader's convenience. 

Consider a bracket [E^/;?] containing the /'th root of some function | (E)[= M\2(E), A(E), etc]. Then, we begin 
each iteration of the root-finding algorithm by computing 

E^(Ef)-Efi}^) 

1 fi*Z(Ef)-&<Z(Ej) ' 1 ' 

where # L and # R are weights necessary to ensure that the regula falsi method converges superlinearly. Then, we 
proceed as follows: 

Tf n^C. n^R. th / If sgnK(£f ° ld )] = sgn[|(£f )]. then # L = I, else # L = l;\ 
If sgn[£ (E) )] = sgn[£ (E) )}, then | £? = ^ R = j h 

m f r/^c,, th / If sgn[§(^' old )] = sgn[^^)], then ,? R = i else t? R =l;l (18) 
Else if sgn[^ = sgn[£ then i >>, 

[ E j ~ E j ' * — L J 

Else { Iteration has converged or bracket was lost. }, 

where £ , ^ old j s the result of Eq. §T7) from the previous iteration. After the conditional statements in Eq. ( fT8| l are 
evaluated, we have a new bracket for the root Ej of %(E), and the algorithm continues by generating a new guess via 
Eq. ( fF7| . The iteration is continued until the bracket is "small enough," i.e., \Ej — Ej\ < e, where e is a tolerance, 
which we took to be 10~ 12 when computing the zeros of A(E) and 10~ 14 when computing the zeros of M21 (E) and 
the ±1 crossings of A(E). 

For the data analysis in the following section, we chose E m j n = — |A| (rj max — fj), where T7 max = max r | T]q(x) \ and 
77 is the (arithmetic) mean of |t}o(jc)|, E max = and Ne = 1000. Note that, while E m ^ n is chosen so that no eigen- 
values exist to the left of it lfT5ll . the value of E max is an empirical estimate based on the spectrum of the data set 
under consideration. Theoretically, E max can be taken to be the Nyquist cutoff, i.e., (n/Ax) 2 , however, the latter is 
very large and impractical. We found that E max = E m [ n + 2\E m i n \ is sufficiently large for most problems, when there is 
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Range (km) x (km) 



Fig. 2. In the left panel, solitary wave packets, generated by the nonhydrostatic Lamb model II II , are shown traveling on the a, = 22.0 isopycnal 
in the (simulated) Yellow Sea. In the right panel, the fourth wave packet from the left (i.e., the range from 159.0 to 164.2 km) is displayed; the zero 
depth in the graph corresponds to a total depth of — hi{= —12 m) on the left. 

no clear choice for the value of E max . Finally, we note that this modified numerical PIST algorithm has been bench- 
marked against the standard examples discussed [ 17], and its results are in exact agreement those of the with previous 
versions/impelementations . 



5. Results and discussion 

In this section, we use the PIST, as described above, to carry out an analysis, along the lines of those in |22 20; 28], of 
a shallow-water internal solitary wave train from the Yellow Sea simulation studies of Chin-Bing et al. [6| and Warn- 
Varnas et al. 11251 . The complete data set and a zoom of the wave packet that we study [i.e., the initial condition r\()(x) 
for the KdV equation] are shown in Fig. [2] Furthermore, the physical constants needed to compute the parameters in 
Eq. Q are taken to be 

g = 9.81m/s 2 , /ii = 12m, h 2 = 58 m, p x = 1020 kg/m 3 , p 2 = 1024.9 kg/m 3 , (19) 

in accordance with those used in 161251 . 

These values result in a nonlinearity-to-dispersion ratio A = —0.00014244, showing that dispersive effects are 
significant, as required for the validity of the KdV-based model. In addition, we can compute an Ursell number (see, 
e.g., EOlO . which is defined as 



16k 2 \ I12 ) \h2, 

for the wavetrain to determine its "degree of nonlinearity." For the wave packet shown in the right panel of Fig. [2] 
we have Ur = 26.8323, where Ur ~ 1 is the typical upper limit of linear theory. Therefore, given the latter values of 
the nonlinearity-to-dispersion ratio and the Ursell number, it is clear that the wave train under consideration here is 
quite nonlinear, with dispersive effects dominating. Also, it is interesting to note that the wave packets in Fig.[2](left) 
are reminiscent of those that emerge from harmonic initial data in the zero-dispersion limit of the KdV equation [9|. 
Indeed, the dimensionless dispersion coefficient is L~ 3 , which is a very small quantity for the internal ocean waves. 

The final test needed to determine the validity of the Osborne-Burch KdV-based model of internal waves for this 
data set is to check whether the fundamental assumption of long waves over shallow water, outlined at the end of 
Section|2j is justified. Clearly, the waves are "long" as we have h\/W ~ 12/200 = 0.06 <C 1, but the "shallowness" 
assumption is barely satisfied as [max x Tjo(x)]//ii ~ 7/12 ss 0.6. In fact, it is easy to see that the latter assumption is 
not valid for the other wave packets featured in the left panel of Fig. |2j they are "too nonlinear" to be described by the 
KdV-based model. This serves to outline the extent to which the nonlinear Fourier analysis is applicable, in general, 
since it is based on the KdV equation, therefore the evolution of the data should be governed by the latter also. 

Nonetheless, even though the wave packet under consideration is barely within the range of validity of the internal 
solitary waves model described in Section[2j we obtain very good results using the scattering transform. To this end, 
in Fig. [3] the Floquet discriminant A(E) is shown. The latter's crossings of the ±1 horizontal lines determine the 
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Fig. 3. The Floquet diagram, which graphically describes the scattering transform spectrum, of the data in Fig. [2] A(£) is plotted logarithmically 
outside of the interval [—1,1], which is denoted by the horizontal dashed lines. The staircase solid line is the accounting function of the Floquet 
discriminant, i.e., S[A(E)], 




0.02 0.03 
2nj/L(m-') 




Fig. 4. In the left panel is the the PIST spectrum of the data in the right panel of Fig. [2] Squares represent the moduli m ; and circles represent the 
amplitudes Aj of the hyperelliptic oscillation modes. In the right panel is the (ordinary) Fourier transform spectrum of the same data set. Note that 
the vertical scale is logarithmic. 



main spectrum eigenvalues {<oj}jLi ■ The dashed vertical line is the spectral reference level E re f, which delineates the 
solitons part of the spectrum from the radiation (and other "less nonlinear" waves) part. Using the spectral reference 
level, we compute the reference depth rj re f = —E Ie f/X = —3.29767 m upon which the solitons propagate. The latter 
is represented in the right panel of Fig.[2]by the dashed horizontal line. Physically, the reference depth of the solitons 
can be understood as follows: in the absence of radiation, rj re f would be the depth at which the unperturbed isopycnal, 
i.e., the the interface between the two layers (recall Fig. |TJ>, is located. 

In Fig.|4](left) the amplitudes and moduli of the nonlinear oscillations modes are plotted versus the (commensurable) 
wavenumbers kj = 2% j/L. For the data set considered here, L — 5203.88 m, M — 342 and N — 36. The dashed vertical 
line corresponds to the wavenumber k le f = 2tzN so ]/L — 0.0132814 m , where iV so i = 11 is the number of solitons. It 
delineates the soliton part of the spectrum from the radiation part, just like the vertical dashed line in Fig. [3] It is clear 
that the spectrum corresponding to the data in the right panel of Fig. |2]consists of eleven solitons. Moreover, there 
is an assortment of moderately nonlinear waves and radiation; however, the amplitudes of the latter modes are very 
small. Therefore, the solitons are the dominant part of the spectrum, and the remaining oscillations do not contribute 
significantly to the physics at hand. 
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Another important finding is that there are "hidden" solitons in the data set. That is to say, in Fig. ^ (right) there 
are clearly four soliton-looking waves, however, the scattering transform spectrum, presented in Fig. B] (left), shows 
there are eleven solitons present. Though the first four are clearly the largest ones, and therefore perceptible to the 
naked eye, there are seven other ones with much smaller amplitudes, and likely larger "wavelengths" (widths), that are 
indiscernible in the data set. This kind of information, which cannot be obtained by means of the (ordinary) Fourier 
transform, is important when studying the effects of internal solitary waves on underwater acoustic signals B24I27I . 

Also in Fig. |4] (right), we have plotted the (ordinary) Fourier spectrum of the data set under consideration. Note 
that we cannot easily read off from the plot the number and type of waves that evolve out of the data, like we did 
above with the scattering transform spectrum. Furthermore, the (ordinary) Fourier spectrum suggest that the most 
energetic modes are in the range, roughly, k 6 [0,0.05]. Meanwhile, the scattering transform spectrum shows that 
the most energetic modes are in the range k £ [0,k le f\ « [0,0.014], which is a fourth of the range predicted by the 
(ordinary) Fourier transfrom. Given the nonlinear nature of the phenomenon the data set describes, we can expect the 
latter wavenumber range to be much more accurate than the former. This ability of the PIST to "pinpoint" the solitons' 
wavenumbers is crucial in studying the acoustic energy transfer (mode conversion) caused by internal waves in the 
ocean 1241271 . 

At this point, it is important to note that the simulated data analyzed in this paper is "clean," in the sense that any 
experimental noise and competing physical effects, other than the generation of solitary solitary waves by the tidal 
flow over bottom topography, are absent. Consequently, this analysis serves as a benchmark for the application of the 
PIST and Osborne's nonlinear Fourier analysis to internal solitary wave data. Now that the limits of applicability have 
been outlined, and the various computational challenges identified and resolved, the PIST is ready for use in the study 
of experimental data, such as the one collected in Strait of Messina and studied by Warn-Varnas et al. Il24l . 

6. Perspective and outlook 

Finally, we mention some possible avenues of future work. An important step towards the wide use of the PIST 
in the oceanographic community is the critical evaluation and benchmarking of the latter against the standard data 
analysis techniques currently in use, i.e., the windowed Fourier and the continuous wavelet transforms. Such work is 
underway by Hawkins et al. ifTOl . Another interesting problem is the development of a "windowed" PIST. That is to 
say, for internal solitary wave trains (such as the one shown in the left panel for Fig. [2] for example) we can compute 
the scattering transform spectrum of each wave packet. Thus, we would have a sequence of spectra illustrating the 
"aging" process the solitons undergo. Much in the same way the windowed Fourier transforms allows for the simul- 
taneous time and frequency resolution of features, a windowed PIST would allow for the identification of constants 
of the motion, i.e., the "nature" of the soliton evolution. The first step in this direction was taken by Zimmerman and 
Haarlemmer l28l . who used this type of "stopwatch" approach to obtain the leading-order invariant of motion of their 
data. Furthermore, such an approach would enable the study of the process by which solitons are created by tidal flow 
over topography, along the lines of the work of Osborne et al. fl9l . 

Another issue that deserves further study is the numerical algorithm for finding the eigenvalues of Eq. Q. The 
method employed herein, like the original automatic algorithm of Osborne fl3], relies on the fact that the data set 
(i.e., the potential of the Schrodinger equation) is piecewise constant to obtain Eq. (13) . This, of course, results in 
a formally first-order method for the spectral eigenvalues. In this regard, there is a promising new approach due to 
Deconinck et al. |8|. The idea is to use a Fourier spectral method in conjunction with Floquet theory to compute the 
eigenvalues (and eigenfunction) of (linear) differential operators such as J^f in Eq. Q. This approach offers greater 
flexibility and improved accuracy, therefore its utility in the nonlinear Fourier analysis of physical wave data should 
be established. 
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